################################################################################
######################20.mlm_single_int.R##############
################################################################################

pols <- readRDS('pols_cwg')

#####Gini:income dif.#####

#install.packages('lme4')
library(lme4)
#install.packages('interplot')
library(interplot)

gg1 <- lmer(scale_econsdw ~ scale_ggini_inc1 +
              scale_socioeconomicsal +
              scale_ggini_inc1 *scale_dispgini +
              scale_dispgini +
              scale_enep +
              scale_galtansdw +
              scale_mag +
              scale_unemp + (1 |  eastwest/ country.name) , 
            data = pols,
            na.action = na.omit,
            REML = T)


pols_rug <- readRDS('pols_bjps') %>%
  dplyr::select(scale_econsdw,scale_dispgini, 
                scale_ggini_inc1,
                scale_socioeconomicsal,
                scale_galtansdw ,
                scale_enep ,
                scale_unemp) %>%
  na.omit

Figure_6 <- interplot(gg1, 'scale_dispgini', 'scale_ggini_inc1', ci = .95) +
  geom_hline(yintercept = 0, color = 'red', linetype = 'dashed')+
  theme_bw() +
  labs(title = '', x = 'Partisan Income Differenetiation', y = 'Marginal Effect of Disposable Income Inequality')+
  geom_rug(data = pols_rug, aes(x = scale_ggini_inc1), inherit.aes = F)

#ggsave('figure_4.pdf', Figure_6, width = 10, height = 4, units = "in")

gg2 <- lmer(scale_econsdw ~ scale_ggini_inc2 +
              scale_socioeconomicsal +
              scale_ggini_inc2 *scale_dispgini +
              scale_dispgini +
              scale_enep +
              scale_galtansdw +
              scale_mag +
              scale_unemp + (1 | eastwest/country.name) , 
            data = pols,
            na.action = na.omit,
            REML = F)




gg3 <- lmer(scale_econsdw ~ scale_ggini_inc3 +
              scale_socioeconomicsal +
              scale_ggini_inc3 *scale_dispgini +
              scale_dispgini +
              scale_enep +
              scale_galtansdw +
              scale_mag +
              scale_unemp + (1 |  eastwest/country.name) , 
            data = pols,
            na.action = na.omit,
            REML = F)




gg4 <- lmer(scale_econsdw ~ scale_ggini_inc4 +
              scale_socioeconomicsal +
              scale_ggini_inc4 *scale_dispgini +
              scale_dispgini +
              scale_enep +
              scale_galtansdw +
              scale_mag +
              scale_unemp + (1 | eastwest/country.name) , 
            data = pols,
            na.action = na.omit,
            REML = F)




gg5 <- lmer(scale_econsdw ~ scale_ggini_inc5 +
              scale_socioeconomicsal +
              scale_ggini_inc5 *scale_dispgini +
              scale_dispgini +
              scale_enep +
              scale_galtansdw +
              scale_mag +
              scale_unemp + (1 |  eastwest/country.name) , 
            data = pols,
            na.action = na.omit,
            REML = F)



#####Point Estimates with gini:ggini#####

#Model Fit Statistics#

r1o <- r.squaredGLMM(gg1)
r2o <- r.squaredGLMM(gg2)
r3o <- r.squaredGLMM(gg3)
r4o <- r.squaredGLMM(gg4)
r5o <- r.squaredGLMM(gg5)

r.o <- rbind(r1o, r2o, r3o, r4o, r5o)
print(mean(r.o[,1]))
print(mean(r.o[,2]))

fit1 <- c(summary(gg1)$AICtab)
fit2 <- c(summary(gg2)$AICtab)
fit3 <- c(summary(gg3)$AICtab)
fit4 <- c(summary(gg4)$AICtab)
fit5 <- c(summary(gg5)$AICtab)

fit.stats <- rbind(fit1, fit2, fit3, fit4, fit5)

print(mean(fit.stats[,1])) #AIC
print(mean(fit.stats[,2])) #BIC
print(mean(fit.stats[,3])) #Log-Likelihood

#SD of Random Intercept from model summaries above#

print(mean(c(0.0000 , 0.0000, 5.474e-11, 0.000, 0.0000)))#Region
print(mean(c(0.8537,0.8193, 8.341e-01 , 0.8352, 0.8183)))#Country

##Regression Table and Dot Plot##

m1o <- as.data.frame(summary(gg1)$coefficients) 

m2o <- as.data.frame(summary(gg2)$coefficients) %>%
  dplyr::rename(Estimate2 = Estimate, 
                `Std. Error2` = `Std. Error`,
                `t value2` = `t value`)
m3o <- as.data.frame(summary(gg3)$coefficients) %>%
  dplyr::rename(Estimate3 = Estimate, 
                `Std. Error3` = `Std. Error`,
                `t value3` = `t value`)
m4o <- as.data.frame(summary(gg4)$coefficients) %>%
  dplyr::rename(Estimate4 = Estimate, 
                `Std. Error4` = `Std. Error`,
                `t value4` = `t value`)
m5o <- as.data.frame(summary(gg5)$coefficients) %>%
  dplyr::rename(Estimate5 = Estimate, 
                `Std. Error5` = `Std. Error`,
                `t value5` = `t value`)

dot.data.o <- cbind(m1o,m2o, m3o,m4o,m5o)

dot.data.o <- dot.data.o[ , order(names(dot.data.o))]

dot.table.o <- dot.data.o %>%
  transmute(pe = rowMeans(dplyr::select(dot.data.o, Estimate, Estimate2, Estimate3, Estimate4, Estimate5)),
            se = rowMeans(dplyr::select(dot.data.o, `Std. Error5`, `Std. Error2`, `Std. Error3`, `Std. Error4`, `Std. Error5`)),
            t.value = rowMeans(dplyr::select(dot.data.o, `t value`, `t value2`, `t value3`, `t value4`, `t value5`))) %>% 
  mutate(lwr = pe - se*1.96) %>%
  mutate(up = pe + se*1.96)

dot.table.o$Predictor <- c("Intercept", "Income.Dif", "Socio-Economic Salience", "Disposable Income Inequality (Gini)", "ENEP", "GALTAN Polarization", "District Magnitude", "Unemployment", "Income.Dif:Gini")

dot.plot.o <- dot.table.o

dot.plot.o$Predictor <- factor(dot.plot.o$Predictor, 
                               levels = c( "Income.Dif", "Socio-Economic Salience", "Disposable Income Inequality (Gini)",  "GALTAN Polarization","ENEP", "District Magnitude", "Unemployment", "Income.Dif:Gini", "Intercept"))

dot.plot.o$Predictor <- factor(dot.plot.o$Predictor, levels=rev(levels(dot.plot.o$Predictor)))

dot.plot.o %>% ggplot(aes(x = Predictor, y = pe)) + 
  geom_point() + 
  geom_hline(aes(yintercept = 0), linetype = 'dashed') +
  geom_linerange(aes(ymin = lwr, ymax = up))+
  coord_flip() +
  labs(title = "Point Estimates of Model Regressing Predictors on Economic Polarization (n = 82)", y = "Point Estimates (95% Confidence Interval)", x = "Predictor")+
  theme_bw()

rownames(dot.table.o) <- dot.table.o$Predictor
dot.table.o <- dot.table.o[,-8]

print(xtable::xtable(dot.table.o, digits = 3))

#####Gini: sal#####
gs1 <- lmer(scale_econsdw ~ scale_ggini_inc1 +
              scale_socioeconomicsal +
              scale_dispgini *scale_socioeconomicsal +
              scale_dispgini +
              scale_enep +
              scale_galtansdw +
              scale_mag +
              scale_unemp+ (1 |  eastwest/country.name) , 
            data = pols,
            na.action = na.omit,
            REML = F)

Figure_7 <- interplot(gs1, 'scale_dispgini', 'scale_socioeconomicsal') +
  geom_hline(yintercept = 0, color = 'red', linetype = 'dashed') +
  theme_bw() +
  labs(title = '',
       x = 'Salience of the Economy', 
       y = '')+
  geom_rug(data = pols_rug, aes(x = scale_socioeconomicsal), inherit.aes = F)

mlm_single_int <- grid.arrange(Figure_6, Figure_7, nrow = 1, top = 'Marginal Effect of Inequality from Mulitlevel Models')

#ggsave("figure_A10.pdf", plot = mlm_single_int, width = 10, height = 4, units = "in")

